global main "C:/Users/silvio/Documents/CVR/reply_to_mvb"
global excel "$main/excel"  
global code "$main/code"  
global data "$main/data"  

global datafile "datacvr" 
global name "ballman_mse_nine"


capture program drop simulmultinomial
program def simulmultinomial
args matrix_in  matrix_out 
matrix mtot=one8*`matrix_in'
matrix p=`matrix_in'/mtot[1,1]
scalar p1=p[1,1]
scalar p2=p[2,1]+p1
scalar p3=p[3,1]+p2
scalar p4=p[4,1]+p3
scalar p5=p[5,1]+p4
scalar p6=p[6,1]+p5
scalar p7=p[7,1]+p6
local ntot=mtot[1,1]
matrix `matrix_out'=J(8,1,0)
forval ijota=1/`ntot' {
scalar isimu =  irecode(runiform(),p1,p2,p3,p4,p5,p6,p7)+1
matrix `matrix_out'[isimu,1]=`matrix_out'[isimu,1]+1
}
end


capture program drop smoothglm
program def smoothglm
args matrix_out 
by j perpe: replace y=max(ntotsimul-y[_N-1]-y[_N-2]-y[_N-3]-y[_N-4]-y[_N-5]-y[_N-6]-y[_N-7],0) if _n==_N &  perpe==perpeh & i==jhat 
glm y u100 u010 u110 u001 u101 u011 u000 if perpe==perpeh & i==jhat , nocons family(poisson) link(log) iterate(20)
qui predict `matrix_out' if perpe==perpeh & i==jhat 
qui replace `matrix_out'=round(`matrix_out')
mkmat `matrix_out' if perpe==perpeh & i==jhat 
drop `matrix_out'
end



matrix drop _all
set more off
set matsize 1392
set seed 98034


cd $code


//User defined macros end

use "${data}/${datafile}" , replace
*n perpe i k j y  (i=j)
by j perpe: replace y=. if _n==_N
by j perpe: replace k=8 if _n==_N

g u100= k==1 | k==3 | k==5 | k==7
g u010= k==2 | k==3 | k==6 | k==7
g u110= k==3 | k==7
g u001= k==4 | k==5 | k==6 | k==7
g u101= k==5 | k==7
g u011= k==6 | k==7
g u111= k==7
g u000=1 

mkmat y


local nsimu=1000
local nstrata=9
local hyponumber=2


forval jj = 1/2 {       // Variate N E  1, N S 2

scalar perpemove=`jj'

if perpemove==1 {
local kk= "State"
}
else if perpemove==2 {
local kk= "Shining_Path"
}

**************************************************************
matrix am=J(8,1,0)
matrix one7=J(1,7,1)
matrix one8=J(1,8,1)

matrix truenobs=J(3,3,0)
**************************************************************



scalar perpemovenot=3-perpemove //Move not
scalar icounter=1
*********************************
//Strata
*********************************
foreach ii in 32 36 25 35 51 14 11 47 48 {  /* begin loop strata */
scalar jhat=`ii'
di jhat

*******************************
//nobs, N direct, indirect
*******************************
scalar perpeh=1
scalar jj1=(jhat-1)*8*3+(perpeh-1)*8+1
scalar jj2=jj1+6
matrix a=y[jj1..jj2,1]              /*original data*/
run "$code/analya"
matrix m_mate=mijk[1..8,modeli]           /*estimated parameters*/
matrix ae=a

scalar perpeh=2
scalar jj1=(jhat-1)*8*3+(perpeh-1)*8+1
scalar jj2=jj1+6
matrix a=y[jj1..jj2,1]              /*original data*/
run "$code/analya"
if modeli>0{
matrix m_mats=mijk[1..8,modeli]           /*estimated parameters*/
}
matrix as=a

*E+S a
matrix a=ae+as                   /*5. Simulated data of sum of E+S */
run "$code/analya"
if modeli>0{
matrix m_mates=mijk[1..8,modeli]           /*estimated parameters*/
}

matrix m_matei=m_mates-m_mats                /*indirect parameters*/
matrix m_matsi=m_mates-m_mate                /*indirect parameters*/

matrix m_mate7=m_mate[1..7,1] 
matrix m_mats7=m_mats[1..7,1] 
matrix m_matei7=m_matei[1..7,1] 
matrix m_matsi7=m_matsi[1..7,1] 

*******************************
//nobs, N direct, indirect end
*******************************

//Direct observations for E,S and ES
//Observations
matrix truenobs[1,1]=one7*m_mate7 
matrix truenobs[2,1]=one7*m_mats7
matrix truenobs[3,1]=truenobs[1,1]+truenobs[2,1]
//Totals Direct
matrix truenobs[1,2]=one8*m_mate 
matrix truenobs[2,2]=one8*m_mats 
matrix truenobs[3,2]=truenobs[1,2]+truenobs[2,2]
//Totals Indirect
matrix truenobs[1,3]=one8*m_matei 
matrix truenobs[2,3]=one8*m_matsi 
matrix truenobs[3,3]=one8*m_mates 

*******************************
//Store True Counts end
*******************************


//Bounds based on Obs, Nin
*******************************
//Move
scalar boundmin=min(truenobs[perpemove,2],truenobs[perpemove,3])
scalar boundmax=max(truenobs[perpemove,2],truenobs[perpemove,3])
scalar  ntruevalue=boundmin

scalar increase=round((boundmax-boundmin)/(`hyponumber'-1))

//Move not
scalar ntruevaluefixed=truenobs[perpemovenot,2]
scalar ntotsimul=ntruevaluefixed
scalar perpeh=perpemovenot

smoothglm m_fix       /*estimated parameters*/
************************************

//Bounds end
*******************************
*********************************
//Loop of Hypothetical values
*********************************
scalar perpeh=perpemove
******************Hypothetical true values for varying
forval itruen = 1/`hyponumber'{

//Varying party
scalar ntotsimul=ntruevalue
************************************

smoothglm m_var      /*estimated parameters. True values for simulation*/

*********************************
//Simulations
*********************************
forval inumi = 1/`nsimu'{              /*inumi=number of simulations*/

matrix nobsxe=J(5,3,0)
matrix nxevari=J(2,1,0)
matrix msevari=J(1,1,0)
matrix msefixvar=J(1,1,0)



**********Simulation********************************************
// Fixed
simulmultinomial m_fix a_fix
matrix a=a_fix[1..7,1]
matrix nobsxe[1,1]=one7*a
matrix nobsxe[2,1]=one8*a_fix  
*********Estimation**************************************************
run "$code/analya"
if modeli>0{
matrix nobsxe[3,1]=round(nmv[modeli, 1])   				/*nctotal*/
matrix nobsxe[4,1]=nmv[modeli, 2] 				     /*vnctotal*/
matrix nobsxe[5,1]=1				     //valid
} /*modeli*/
scalar modeli_fix=modeli
***************************************************************************

**********Simulation********************************************
// Varying

simulmultinomial m_var a_var

matrix a=a_var[1..7,1]
matrix nobsxe[1,2]=one7*a
matrix nobsxe[2,2]=one8*a_var  
*********Estimation**************************************************
run "$code/analya"
if modeli>0{
matrix nobsxe[3,2]=round(nmv[modeli, 1])   				/*nctotal*/
matrix nobsxe[4,2]=nmv[modeli, 2] 				     /*vnctotal*/
matrix nobsxe[5,2]=1				     //valid
} /*modeli*/
scalar modeli_var=modeli
***************************************************************************

*Fixed+Varying
matrix am=a_fix+a_var                   /*5. Simulated data of sum of E+S */
matrix a=am[1..7,1]
matrix nobsxe[1,3]=one7*a
matrix nobsxe[2,3]=one8*am  
**********************************matrix xe[`inumi',1]=one7*am  /*1. Simulated for Fixed */
*Estimate on simulated data
run "$code/analya"
if modeli>0{
matrix nobsxe[3,3]=round(nmv[modeli, 1])   				/*nctotal*/
matrix nobsxe[4,3]=nmv[modeli, 2] 				     /*vnctotal*/
matrix nobsxe[5,3]=1				     //valid
} /*modeli*/
***************************************************************************

//Varying
if modeli>0 &  modeli_fix>0 {                   //Estimates have to be greater than observed data
if nobsxe[3,3]- nobsxe[2,1]-nobsxe[1,2]>0{               //Estimates have to be greater than observed data
matrix nxevari[2,1]=1
}
}

//mse var 
matrix msefixvar[1,1]= (nobsxe[3,2]- nobsxe[2,2])*(nobsxe[3,2]- nobsxe[2,2])* nobsxe[5,2]

matrix nxevari[1,1]= (nobsxe[3,3]- nobsxe[3,1])*nxevari[2,1]
matrix msevari[1,1]= (nobsxe[3,3]- nobsxe[3,1]-nobsxe[2,2])*(nobsxe[3,3]- nobsxe[3,1]-nobsxe[2,2])*nxevari[2,1]



****************************************************************
if `inumi'==1{
matrix nobsxetot=nobsxe
matrix nxevaritot=nxevari
matrix msevaritot=msevari
matrix msefixvartot=msefixvar

}
else{
matrix nobsxetot=nobsxetot+nobsxe
matrix nxevaritot=nxevaritot+nxevari
matrix msevaritot=msevaritot+msevari
matrix msefixvartot=msefixvartot+msefixvar
}

} /*simulation*/
*********************************
//Simulations end
*********************************

preserve
clear

forval i = 1/3{
forval j = 1/4{
matrix nobsxetot[`j',`i']=nobsxetot[`j',`i']/nobsxetot[5,`i']
}
}

matrix msefixvartot[1,1]=msefixvartot[1,1]/nobsxetot[5,1]

matrix nxevaritot[1,1]=nxevaritot[1,1]/nxevaritot[2,1]
matrix msevaritot[1,1]=msevaritot[1,1]/nxevaritot[2,1]


di `inumi',`itruen',icounter,jhat

matrix res= jhat,ntruevalue
matrix res_nvar=nobsxetot[3,2]
matrix res_msevar=msefixvartot[1,1]/1000

matrix res_nvari=nxevaritot[1,1]
matrix res_msevari=msevaritot[1,1]/1000

if icounter==1{
matrix res_tot_var=res,res_nvar,res_msevar,res_nvari,res_msevari
}
else {
matrix res_tot_var=res_tot_var\res,res_nvar,res_msevar,res_nvari,res_msevari
}

scalar icounter=icounter+1			  
scalar ntruevalue=ntruevalue+increase

restore
} /* end loop true N hypothetical true values*/
*********************************
//Hypothetical values end
*********************************

} /* end loop strata*/
*********************************
//Strata end
*********************************

*matrix list res
preserve
clear

*************************************************************
*Export output
*************************************************************


svmat res_tot_var, names(v)

ren v1 stratum		
ren v2 ntruevalue

ren v3 nobstot_direct
ren v4 mse_direct

ren v5 nobstot_indirect
ren v6 mse_indirect


export excel * ///
using "$excel/${name}.xls", sheet("results_`kk'") sheetmodify firstrow(varlabels) missing(".")
clear


restore
} //j

clear

